###########################################
### MAIN                                ###
###########################################

# Data preparation --------------------------------------------------------

# load packages

pckgs <-
  c(
    "raster",
    "grid",
    "gridExtra",
    "dplyr",
    "ggplot2",
    "devtools",
    "Hmisc",
    "glmmTMB",
    "sjmisc",
    "sjlabelled",
    "sjPlot",
    "sjstats",
    "ggeffects",
    "parameters",
    "scales")

lapply(pckgs, require, character.only = TRUE)

# set wd

setwd("C:/Users/majacob/polybox/Statistics/R/2021_party_income_states/replication")

# loading data

data_maps <-
  read.csv(file = "Germany_finance_states_inc_heat_map_2009-2017.csv", na.strings =
             c("", "NA"))

data_regression  <-
  read.csv(file = "Germany_finance_states_inc_analysis_2009-2017.csv", na.strings =
             c("", "NA"))

# Figure 1 ----------------------------------------------------------------

# Compute shares

data_regression$total_income <- data_regression$i_fees+data_regression$i_man+data_regression$i_don_ind+data_regression$i_don_cor+data_regression$i_sub
data_regression$share <- (data_regression$i_don_cor+data_regression$i_don_ind)/data_regression$total_income
data_regression$share_individual <- (data_regression$i_don_ind)/data_regression$total_income
data_regression$share_corporation <- (data_regression$i_don_cor)/data_regression$total_income
data_regression$total_income_in_100k <- data_regression$total_income/100000

# Rename states

data_regression$state <- memisc::recode(as.factor(data_regression$state),"Brandenburg" <- "BB","Berlin"<- "BE",
                             "Baden-Württemberg" <- "BW",
                             "Bavaria" <- "BY",
                             "Bremen" <- "HB",
                             "Hesse" <- "HE",
                             "Hamburg" <- "HH",
                             "Meckl.-Vorpommern" <- "MV",
                             "Lower Saxony" <- "NI",
                             "North Rhine-Westph." <- "NW",
                             "Rhineland-Palatinate" <- "RP",
                             "Schleswig-Holstein" <- "SH",
                             "Saarland" <- "SL",
                             "Saxony" <- "SN",
                             "Saxony-Anhalt" <- "ST",
                             "Thuringia" <- "TH")

data_regression$Party <- data_regression$party

# Build graphs

share_total <- ggplot(data = data_regression ,
                      aes(x = year,
                          y = share,
                          shape=Party)) + theme_bw() + xlab("Year") + geom_point() + ylab("Share of Donations in Total Income") + facet_wrap( ~ state, nrow = 3)  + scale_x_continuous(limits = c(2009, 2017), breaks = c(2009,2011,2013,2015,2017)) + geom_line()


share_total
share_individual <- ggplot(data = data_regression ,
                           aes(x = year,
                               y = share_individual,
                               shape=Party)) + theme_bw() + xlab("Year") + geom_point() + ylab("Share of Individual Donations in Total Income") + facet_wrap( ~ state, nrow = 3)  + scale_x_continuous(limits = c(2009, 2017), breaks = c(2009,2011,2013,2015,2017)) + geom_line()


share_corporation <- ggplot(data = data_regression ,
                            aes(x = year,
                                y = share_corporation,
                                shape=Party)) + theme_bw() + xlab("Year") + geom_point() + ylab("Share of Corporative Donations in Total Income") + facet_wrap( ~ state, nrow = 3)  + scale_x_continuous(limits = c(2009, 2017), breaks = c(2009,2011,2013,2015,2017)) + geom_line()


absolute_income <- ggplot(data = data_regression ,
                          aes(x = year,
                              y = total_income_in_100k,
                              shape=Party)) + theme_bw() + xlab("Year") + geom_point() + ylab("Total Income in 100,000k Euro") + facet_wrap( ~ state, nrow = 3,scales="free_y")  + scale_x_continuous(limits = c(2009, 2017), breaks = c(2009,2011,2013,2015,2017)) + geom_line()

# Plot

pdf(file = "Figure_1.pdf",
    width = 10,
    height = 6,
    onefile=F)

share_total + theme(axis.text.x = element_text(size = 6, angle = 45,hjust = 1.25, vjust = 1.5))

dev.off()

# Figure 2 ----------------------------------------------------------------

cdu <- make_map_individual_donations(data=data_maps[which(data_maps$party == 'CDU'), ],title="CDU")
spd <- make_map_individual_donations(data=data_maps[which(data_maps$party == 'SPD'), ],title="SPD")
fdp <- make_map_individual_donations(data=data_maps[which(data_maps$party == 'FDP'), ],title="FDP")
left <- make_map_individual_donations(data=data_maps[which(data_maps$party == 'Left'), ],title="Left")
greens <- make_map_individual_donations(data=data_maps[which(data_maps$party == 'Greens'), ],title="Greens")
afd <- make_map_individual_donations(data=data_maps[which(data_maps$party == 'AfD'), ],title="AfD")

# Plot

pdf(file = "Figure_2.pdf",
    width = 10,
    height = 6,
    onefile=F)

raster::update(c(cdu, spd, fdp, left, greens,afd), layout = c(3,2), as.table = TRUE)

dev.off()

# Figure 3 ----------------------------------------------------------------

cdu <- make_map_corporate_donations(data=data_maps[which(data_maps$party == 'CDU'), ],title="CDU")
spd <- make_map_corporate_donations(data=data_maps[which(data_maps$party == 'SPD'), ],title="SPD")
fdp <- make_map_corporate_donations(data=data_maps[which(data_maps$party == 'FDP'), ],title="FDP")
left <- make_map_corporate_donations(data=data_maps[which(data_maps$party == 'Left'), ],title="Left")
greens <- make_map_corporate_donations(data=data_maps[which(data_maps$party == 'Greens'), ],title="Greens")
afd <- make_map_corporate_donations(data=data_maps[which(data_maps$party == 'AfD'), ],title="AfD")

# Plot

pdf(file = "Figure_3.pdf",
    width = 10,
    height = 6,
    onefile=F)

raster::update(c(cdu, spd, fdp, left, greens,afd), layout = c(3,2), as.table = TRUE)

dev.off()

# Table 1 -----------------------------------------------------------------

data <- data_regression

# prepare data

data$i_don <- data$i_don_ind + data$i_don_cor # building overall donation variable
data$inh_1000 <- data$inh * 1000 # getting overall number of inhabitants
data$GDP_1000000 <- data$gdp * 1000000 # getting overall amount of GDP

# set variables to factors

data$fed_opp <- car::recode(data$federal_gov, "0=1;1=0", as.factor = TRUE)
data$state_gov <- as.factor(data$state_gov)
data$state_el <- as.factor(data$state_el)
data$federal_el <- as.factor(data$federal_el)

# round DVs (for negative binomial model)

data$i_don <- round(data$i_don, digits = 0)
data$i_don_ind <- round(data$i_don_ind, digits = 0)
data$i_don_cor <- round(data$i_don_cor, digits = 0)

# rescale IVs

data$inh_1000sd <-  scales::rescale(data$inh_1000,to = c(0,100))
data$inh_1000sd <- round(data$inh_1000sd, digits = 0)

data$GDP_1000000sd <-  scales::rescale(data$GDP_1000000,to = c(0,100))
data$GDP_1000000sd <- round(data$GDP_1000000sd, digits = 0)

# transform ENP

data$ENPx10 <- data$ENP * 10 # multiply ENP times 10 to receive integers
data$ENPx10 <- round(data$ENPx10, digits = 0)

# build dataframe without AfD observations

datamAfD <- data[!(data$party == "AfD"), ]

# Model l to 3

model_don_1_1 <-
  glmmTMB(
    i_don ~ ENPx10 + state_gov + fed_opp + state_el + federal_el
    + GDP_1000000sd + inh_1000sd + (1 | state) + (1 | party/state),
    family = nbinom2(link = "log"),
    REML = FALSE,
    data = data
  )

summary(model_don_1_1)

model_don_1_2 <-
  glmmTMB(
    i_don ~ ENPx10 + state_gov + fed_opp + state_el + federal_el + ENPx10 *
      state_el
    + GDP_1000000sd + inh_1000sd + (1 |
                                      state) + (1 |
                                                  party/state),
    family = nbinom2(link = "log"),
    REML = FALSE,
    data = data
  )

model_don_2 <-
  glmmTMB(
    i_don ~ ENPx10 + state_gov + fed_opp + state_el + federal_el + ENPx10 *
      state_el
    + GDP_1000000sd  + inh_1000sd + (1 |
                                      state) + (1 |
                                                  party/state),
    family = nbinom2(link = "log"),
    REML = FALSE,
    data = datamAfD
  )

summary(model_don_1_2)

# Model 4

model_ind <-
  glmmTMB(
    i_don_ind ~ ENPx10 + state_gov + fed_opp + state_el + federal_el + ENPx10 *
      state_el
    + GDP_1000000sd + inh_1000sd + (1 |
                                      state) + (1 |
                                                  party/state),
    family = nbinom2(link = "log"),
    REML = FALSE,
    data = data
  )

summary(model_ind)


# Model 5

model_cor <-
  glmmTMB(
    i_don_cor ~ ENPx10 + state_gov + fed_opp + state_el + federal_el + ENPx10 *
      state_el
    + GDP_1000000sd + inh_1000sd + (1 |
                                      state) + (1 |
                                                  party/state),
    family = nbinom2(link = "log"),
    REML = FALSE,
    data = data
  )

# Table

pl <- c(
  `(Intercept)` = "Intercept",
  ENPx10 = "ENEPx10",
  state_gov1 = "State government",
  fed_opp1 = "Federal opposition",
  state_el1 = "State election",
  federal_el1 = "Federal election",
  GDP_1000000sd = "State GDP",
  inh_1000sd = "State inhabiants",
  `ENPx10:state_el1` = "(ENEPx10)xState election"
)

tab_model(
  model_don_1_1,
  model_don_1_2,
  model_don_2,
  model_ind,
  model_cor,
  show.ci = FALSE,
  show.p = FALSE,
  transform = NULL,
  show.std = "std2",
  show.zeroinf = TRUE,
  show.est = FALSE,
  show.se = TRUE,
  collapse.se = TRUE,
  digits = 3,
  p.style = "stars",
  dv.labels = c(
    "Model 1 (all donations)",
    "Model 2 (all donations)",
    "Model 3 (all donations)",
    "Model 4 (individual donations)",
    "Model 5 (corporative donations)"
  ),
  pred.labels = pl,
  show.re.var = FALSE,
  show.r2	= FALSE,
  show.icc = FALSE,
  show.loglik	= TRUE,
  show.dev = TRUE
  
)

# Figure 4 ----------------------------------------------------------------

# Predicted counts and building plots

margeff_model_ind <-
  ggeffects::ggpredict(model_ind, c("ENPx10 [all]", "state_el"))
options(scipen = 10000)
margeff_model_ind_plot <-
  plot(margeff_model_ind, colors = "bw", show.legend = TRUE,add.data=T,ci.style="ribbon",jitter=0.3,use.theme=T) + labs(
    x = "ENEP*10",
    y = "Individual donations in euro",
    title = "a)",
    linetype  = "State election"
  ) + scale_linetype_manual(values = 15:17, labels = c("no","yes")) + ylim(0,1000000)

margeff_model_ind_plot

margeff_model_cor <-
  ggeffects::ggpredict(model_cor, c("ENPx10 [all]", "state_el"))
margeff_model_cor_plot <-
  plot(margeff_model_cor, colors = "bw", show.legend = TRUE,add.data=T,ci.style="ribbon",jitter=0.3,use.theme=T) + labs(
    x = "ENEP*10",
    y = "Corporative donations in euro",
    title = "b)",
    linetype  = "State election"
  ) + scale_linetype_manual(values = 15:17, labels = c("no","yes")) + ylim(0,340000)

margeff_model_cor_plot

# Plot

pdf(file = "Figure_4.pdf",
    width = 10,
    height = 5,
    onefile=F)

grid.arrange(margeff_model_ind_plot,margeff_model_cor_plot, nrow=1)

dev.off()